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Abstract Development of the Aarhus adiabatic pulsation 
code started around 1978. Although the main features have 
been stable for more than a decade, development of the code 
is continuing, concerning numerical properties and output. 
The code has been provided as a generally available pack- 
age and has seen substantial use at a number of installations. 
Further development of the package, including bringing the 
documentation closer to being up to date, is planned as part 
of the HELAS Coordination Action. 
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1 Introduction 

The goal of the development of the code was to have a sim- 
ple and efficient tool for the computation of adiabatic os- 
cillation frequencies and eigenfunctions for general stellar 
models, emphasizing also the accuracy of the results. Not 
surprisingly, given the long development period, the sim- 
plicity is now less evident. However, the code offers con- 
siderable flexibility in the choice of integration method as 
well as ability to determine all frequencies of a given model, 
in a given range of degree and frequency. 

The choice of variables describing the equilibrium model 
and os cillations was to a large exte nt inspired bv Dziem- 
bowski lll97ll) . As discussed in Section 12. 11 the equilibrium 
model is defined in terms of a minimal set of dimensionless 
variables, as well as by mass and radius of the model. 
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Fairly extensive documentation of the code, on which 
the present paper i n part is based, is provided with the dis- 
tributi on packageQ. IChristensen-Dalsgaard and Berthomieul 
(Il99ll) provided an extensive review of adiabatic stellar os- 
cillations, emphasizing applications to helioseismology, and 
discusse d many aspects and tests of the A a rhus p ackage, 
whereas IChristensen-Dalsgaard and MullanI d 1994 1) carried 
out careful tests and comparisons of results on polytropic 
models; this includes extensive tables of frequencies which 
can be used for comparison with other codes. 



2 Equations and numerical scheme 

2.1 Equilibrium model 

The equilibrium model is defined in terms of the following 
dimensionless variables: 
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Fi dlnr F^pr 



1 dln/j dlnp 



Fi dlnr dlnr ' 



(1) 



Here r is distance to the centre, m is the mass interior to 
r, R is the photospheric radius of the model and M is its 
mass; also, G is the gravitational constant, p is pressure, p 
is density, and Fi = {dlnp/dlnp),^^^, the derivative being at 
constant specific entropy. In addition, the model file defines 
M and R, as well as central pressure and density, in dimen- 
sional units, and scaled second derivatives of p and p at the 
centre (required from the expansions in the central boundary 

' The package is available at 

|http: // astro . p hys . au . dk/ ^ j cd/ adipack . ii| 
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condition); finally, for models with vanishing surface pres- 
sure, assuming a polytropic relation between p and p in the 
near-surface region, the polytropic index is specified. 

The following relations between the variables defined 
here and more "physical" variables are often useful: 

f2 



P = 



P = 



GM^ x^AfAs dp 
AnR^ A9A3 ' dr 



4nR^ 



xAiA 



1^5 



1^5 



(2) 



We may also express the characteristic frequencies for adi- 
abatic oscillations in terms of these variables. Thus if is 
the buoyancy frequency. Si is the Lamb frequency at degree 
/ and (Oa is the acoustical cut-off frequency for an isothermal 
atmosphere, we have 



Sf = 



GM 



GM 



-A, A 



1/14 



GM .2 _ GMl{l+l)Ai 



GM 



1 GM 



(3) 



(4) 



(5) 



where c is the adiabatic sound speed, and Hp=p/ (gp) is the 
pressure scale height, g being the gravitational acceleration. 
Finally it may be noted that the squared sound speed is given 
by 



2 GM ^2 
c — c 

R 



GM 2M 



(6) 



These equations also define the dimensionless characteristic 
frequencies N, §1 and (Ua as well as the dimensionless sound 
speed c, which are often useful. 



2.2 Formulation of the equations 

As is well known the displacement vector of nonradial 
(spheroidal) modes can be written in terms of polar coor- 
dinates (r, 0,0) as 



5r = Re 



(7) 



1 dY" 



sine 



exp(— ift)f) 



Here y"'(0,(/)) = c/„,/'™(cos 0)exp(im(/)) is a spherical har- 
monic of degree / and azimuthal order m, 9 being co-latitude 
and (/) longitude; P'"{x) is an associated Legendre function, 
and cjin is a suitable normalization constant. Also, a,-, ag, and 
Hip are unit vectors in the r, 0, and directions. Finally, f is 
time and ft) is the angular frequency of the mode. Similarly, 
e.g., the Eulerian perturbation to pressure may be writterQ 



/(r, ,(/), = Re [/ (r)y/" ( ,(/)) exp( -ift)0] 



(8) 



^ I do not here distinguish between the full perturbation and the ra- 
dial amphtude function. 



As the oscillations are adiabatic (and only conservative 
boundary conditions are considered) ft) is real, and the am- 
plitude functions ^h('')j /''('')> etc. can be chosen to be 
real. 

The equations of adiabatic stellar oscillations, in the non- 
radial case, are expressed in terms of the following vari- 
ables0 



yi 



R ' 



VP / ft)^r^ R 
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gr 

- x^— (—\ 
60c \ X J 



(9) 



Here is the perturbation to the gravitational potential. 
Also, we introduce the dimensionless frequency a by 



2 GM 2 



(10) 



corresponding to Eqs|3]-|5] These quantities satisfy the fol- 
lowing equations: 



dx 



(11) 



= [/a+l)-77A]yi + (A-l)3;2 + 7?AB, (12) 

dy3 
ax 



x^ = -AUyi-U^y2 
dx 77 



(13) 
(14) 



-[/(/+ 1) + [/(A - 2) + UVg]y3 + 2(1 - U)y4 



Here 77 =l{l + l)g/io}'^r) = 1(1 + 1 )A| / g^, and the not ation 
is otherwise as defined in Eq. [T] In the iCowlingI d 19411) ap- 
proximation, where the perturbation to the gravitational po- 
tential is neglected, the terms in yi, are neglected in Eqs [TT] 
and[l2]and Eqs [13] and [14] are not used. 

The dependent variables y, in the nonradial case have 
been chosen in such a way that for / > they all vary as 
x'^' for jc — > 0. For large I a considerable (and fundamen- 
tally unnecessary) computational effort would be needed to 
represent this variation sufficiently accurately with, e.g., a 
finite difference technique, if these variables were to be used 
in the numerical integration. Instead I introduce a new set of 
dependent variables by 



--x-'+'y,, 



i = 1,2,3,4 



(15) 



These variables are then 0(1) inx near the centre. They are 
used in the region where the variation in the j, is dominated 

^ The somewhat pecuhar choice of y^, y4 resuhs from the earUer 
use of an unconventional sign convention for <P'; now, as usual, <P' is 
defined such that the perturbed Poisson equation has the form V^<P' = 
AnGp' , where p' is the Eulerian density perturbation. 
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by the behaviour, for x < jCgv, say, where Xgv is deter- 
mined on the basis of the asymptotic properties of the solu- 
tion. This transformation permits calculating modes of arbi- 
trarily high degree in a complete model. 

For radial oscillations only yi and are used, where yi 
is defined as above, and 



The code also has the option of considering truncated 
(e.g., envelope) models although at the moment only in the 
Cowling approximation or for radial oscillations. In this case 
the innermost boundary condition is typically the vanish- 
ing of the radial displacement although other options are 
available. 



Here the equations become 

dy2 _ 
dx 



-yi 



yi +Ay2 



(16) 



(17) 



The equations are solved on the interval [xi,Xs] in x. 
Here, in the most common case involving a complete stellar 
model x\ = e, where £ is a suitably small number such that 
the series expansion around x = is sufficiently accurate; 
however, the code can also deal with envelope models with 
arbitrary xi, typically imposing = at the bottom of the 
envelope. The outermost point is defined by Xs = R^/R where 
/?s is the surface radius, including the atmosphere; thus, typ- 
ically, Xs > 1 . 



2.3 Boundary conditions 

The centre of the star, r = 0, is obviously a singular point of 
the equations. As discusse d, e.g.. by Christensen-Dalsgaard 
et al. (il974J) boundary conditions at this point are obtained 
from a series expansion, in the present code to second sig- 
nificant order. In the general case this defines two conditions 
at the innermost non-zero point in the model. For radial os- 
cillations, or nonradial oscillations in the Cowling approxi- 
mation, one condition is obtained. The surface in a realistic 
model is typically defined at a suitable point in the stellar at- 
mosphere, with non-zero pressure and density. Here the sim- 
ple condition of vanishing Lagrangian pressure perturbation 
is implemented and sometimes used. However, more com- 
monly a condition between pressure perturbation and dis- 
placement is established by matching continuously to the so- 
lution in an isothermal atmosphere extending continuously 
from the uppermost poin t in the m odel Q A very similar con- 
dition was presented by IUnno"e t al. (1989). In addition, in 
the full nonradial case a condition is obtained from the con- 
tinuous match of <P' and its derivative to the vacuum solution 
outside the star. 

In full polytropic models, or other models with vanish- 
ing surface pressure, the surface is also a singular point. 
In this case a boundary condition at the outermost non- 
singular point is obtained from a series expansion, assum- 
ing a near-surface polytrop ic behaviour (see Christensen- 
Dalsgaard and Mullan. 11994 for details). 

^ Note that since the frequency, and other variables, are taken to be 
real this can only be applied for frequencies below the acoustical cut- 
off frequency in the isothermal extension. 



2.4 Numerical scheme 

The numerical problem can be formulated generally as that 
of solving 



dyi V 



(18) = X1«'7WX/W ' for/= 1,...,/, 



(19) 



with the boundary conditions 



£Z?yy,(xi)=0, for /=1,... ,7/2, (20) 



7=1 



for i = 1 



.,1/2. 



(21) 



Here the order I of the system is 4 for the full nonradial 
case, and 2 for radial oscillations or nonradial oscillations in 
the Cowling approximation. This system only allows non- 
trivial solutions for selected values of which is thus an 
eigenvalue of the problem. 

The programme permits solving these equations with 
two basically different techniques, each with some variants. 
The first is a shooting method, where solutions satisfying the 
boundary conditions are integrated separately from the inner 
and outer boundary, and the eigenvalue is found by match- 
ing these solutions at a suitable inner fitting point Xf. The 
second technique is to solve the equations together with a 
normalization condition and all boundary conditions using 
a relaxation technique; the eigenvalue is then found by re- 
quiring continuity of one of the eigenfunctions at an interior 
matching point. 

For simplicity I do not distinguish between y, and y, (cf. 
Section |2.2| l in this section. It is implicitly understood that 
the dependent variable (which is denoted y,) is y, for x < Xev 
and y, for x > Xev The numerical treatment of the transition 
between y, and y, has required a little care in the coding. 



2.5 The shooting method 

It is convenient here to distinguish between 1 = 2 and / = 
4. For 7 = 2 the differential Eqs [19] have a unique (apart 

from normalization) solution yf^ satisfying the inner bound- 
ary conditions [20] and a unique solution y|°' satisfying the 
outer boundary conditions |2T] These are obtained by numer- 
ical integration of the equations. The final solution can then 

be represented as yj = C^'^y^' = C^°^y^°\ The eigenvalue is 
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obtained by requiring that the solutions agree at a suitable 
matching point x = Xf, say. Thus 

dSi"(^f) = c(°)/°'(xf), 

C»j«(xf)=d°)y(°'(xf). (22) 

These equations clearly have a non- trivial solution (CW,C(°)) 
only when their determinant vanishes, i.e., when 



The programme allows the use of two different tech- 
niques for solving the differential equations. One is the stan- 
dard second-order centred difference technique, where the 
differential equations are replaced by the difference equa- 
tions 



yi 



v«+l 



,n+l 



i = l,...,/. (29) 



0. 



(23) 



Equation |23] is therefore the eigenvalue equation. 

For 7 = 4 there are two linearly independent solutions 
satisfying the inner boundary conditions, and two linearly 
independent solutions satisfying the outer boundary condi- 
tions. The former set {jf '^^yf ^'} is chosen by setting 



Here I have introduced a mesh xi = < < ■ ■ • < = Xs 
in X, where A' is the total number of mesh points; = y, (x"), 
and = aji{x"). These equations allow the solution at x = 

x"+' to be determined from the solution at x = x". 

The se cond technique was proposed by Gabriel and 
Noels (197^ here on each mesh interval (x' 
sider the equations 



we con- 



yP(^i) = l: 



dy] 



(«) 



(24) dx 



L-n (n 



for i — 1 



(30) 



and the latter set {y\°'^\y\°'^^ } is chosen by setting 



/r'(^s) = i, 



0, 



^3 (^'s) = l 



(25) 



The inner and outer boundary conditions are such that, given 
yi and ^3, y2 and ^4 may be calculated from them; thus 
Eqs [24] and |25] completely specify the solutions, which are 
obtained by integrating from the inner or outer boundary. 
The final solution can then be represented as 



yj = d''''y 



(M) 



C 



(i,2)^(i,2) 



c 



{o,2)-y(o-2) 



(26) 



with constant coefficients, where a"j = l/2{a"j + a"^^). 
These equations may be solved analytically on the mesh in- 
tervals, and the complete solution is obtained by continuous 
matching at the mesh points. This technique clearly permits 
the computation of solutions varying arbitrarily rapidly, i.e., 
the calculation of modes of arbitrarily high order. On the 
other hand solving Eqs [30] involves finding the eigenvalues 
and eigenvectors of the coefficient matrix, and therefore be- 
comes very complex and time consuming for higher-order 
systems. Thus in practice it has only been implemented for 
systems of order 2, i.e., radial oscillations or nonradial os- 
cillations in the Cowling approximation. 



At the fitting point Xf continuity of the solution requires that 

C('''));f (xf ) + C'^^yj^ (xf) = (27) 2.6 The relaxation technique 

d°-'V;'''\xt)+d°^^V;''\xf) 7 = 1,2,3,4. 
This set of equations only has a non-trivial solution if 



:det< 



(M) (i,2) {0,1) 

yis yii y\,i 

(ia) (i,2) {0,1) 

yis yii y2,i 

ill) (12) {0,1) 

y3s y?,i y3,i 

(ia) (i,2) {0,1) 

y4.f 3'4.f y4,i 



(0-2) 
3'l.f 

Jo,2) 

yi.i 
(0,2) 

(0,2) 

3'4.f 



(28) 



where, e.g., y^^'p 



y^j'^^ (xf). Thus Eq. 



is the eigenvalue 
5] or Eq. |28] is a 



equation in this case. 

Clearly A as defined in either Eq. 
smooth function of a^, and the eigenfrequencies are found 
as the zeros of this function. This is done in the programme 
using a standard secant technique. However, the programme 
also has the option for scanning through a given interval in 
CT^ to look for changes of sign of A, possibly iterating for the 
eigenfrequency at each change of sign. Thus it is possible to 
search a given region of the spectrum completely automati- 
cally. 



If one of the boundary conditions is dropped, the difference 
equations, with the remaining boundary condition and a nor- 
malization condition, constitute a set of linear equations for 
the {y'j} which can be solved for any value of a; this set may 
be solved e fficiently by forwa rd elimination and backsubsti- 
tution (e.g. . .Baker et al.lfT97 ih . with a pr ocedure very similar 
to the so-called Henyey technique (e.g.. lHenvey et al.L[l959l 
see also Christensen-Dalsgaard 2007) used in stellar mod- 
elling. The eigenvalue is then found by requiring that the re- 
maining boundary condition, which effectively takes the role 
of A(a), be satisfied. However, as both boundaries, at least 
in a complete model, are either singular or very nearly sin- 
gular, the removal of one of the boundary conditions tends 
to produce solutions that are somewhat ill-behaved, in par- 
ticular for modes of high degree. This in turn is reflected in 
the behaviour of Zi as a function of a. 

This problem is avoided in a variant of the relaxation 
technique where the difference equations are solved sepa- 
rately for X < Xf and x > Xf, by introducing a double point 
XjT = x"f = x"f+^ = x^ in the mesh. The solution is further- 
more required to satisfy the boundary conditions |20l and 1271 
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a suitable normalization condition (e.g. yi (x^) = 1), and con- 
tinuity of all but one of the variables atx = Xf, e.g., 

yiixf) =yi{x+) , 

y4{x:[) =y4ix^) , (31) 

(when 1 = 2 clearly only the first continuity condition is 
used) We then set 

A=y2iXf)~y2{x+), (32) 

and the eigenvalues are found as the zeros of A , regarded as 
a function of a^. With this definition, A may have singular- 
ities with discontinuous sign changes that are not associated 
with an eigenvalue, and hence a little care is required in the 
search for eigenvalues. However, close to an eigenvalue A 
is generally well-behaved, and the secant iteration may be 
used without problems. 

As implemented here the shooting technique is consid- 
erably faster than the relaxation technique, and so should 
be used whenever possible (notice that both techniques may 
use the difference Eqs|29]and so they are numerically equiv- 
alent, in regions of the spectrum where they both work). For 
second-order systems the shooting technique can probably 
always be used; the integrations of the inner and outer so- 
lutions should cause no problems, and the matching deter- 
minant in Eq.|23]is well-behaved. Vov fourth-order systems, 
however, this needs not be the case. For modes where the 
perturbation to the gravitational potential has little effect on 

the solution, the two solutions y^''^' and y^''^', and similarly 

the two solutions y'^" '' and y''°'^\ are almost linearly depen- 
dent, and so the matching determinant nearly vanishes for 
any value of a^. This is therefore the situation where the 
relaxation technique may be used with advantage. This ap- 
plies, in particular, to the calculation of modes of moderate 
and high degree which are essential to helioseismology. 

2.7 Improving the frequency precision 

To make full use of the increasingly accurate observed fre- 
quencies the computed frequencies should clearly at the very 
least match the observational accuracy, for a given model. 
Only in this way do the frequencies provide a faithful rep- 
resentation of the properties of the model, in comparisons 
with the observations. However, since the numerical errors 
in the computed frequencies are typically highly systematic, 
they may affect the asteroseismic inferences even if they 
are smaller than the random errors in the observations, and 
hence more stringent requirements should be imposed on the 
computations. Also, the fact that solar-like oscillations, and 
several other types of asteroseismically interesting modes, 
tend to be of high radial order complicates reaching the re- 
quired precision. 

The numerical techniques discussed so far are generally 
of second order. This yields insufficient precision in the eval- 
uation of the eigenfrequencies, unless a very dense mesh is 



used in the computation (see also iMoya et all l2007h . The 
code may apply two tec hniques to improve the precisi on. 

One technique (cf. IChristensen -Dalsgaardl Il982h uses 
the fact that the frequency approxim ately satisfies a varia- 
tional principle (' Chandrasekhari fl964l) Fl The variational ex- 
pression may formally be written as 

a^ = al,^E[^f = ^, (33) 

where and are integrals over the equilibrium model 
depending on the eigenfunction, here represented by ^. 
The variational property implies that any error 5S, in ^ in- 
duces an error in that is Thus by substitut- 
ing the computed eigenfunction into the variational expres- 
sion a more precise determination of should resul t. This 
has indeed been confirmed (fch ristensen- DalsgaMl Il982t 
IChristensen-Dalsgaard and Berthomieul 119911: Christensen- 
Dalsgaard and Mullan. ll994l) . 

The second technique uses explicitly that the difference 
scheme |29] which is used by one version of the shooting 
technique, and the relaxation technique, is of second order. 
Consequently the truncation errors in the eigenfrequency 
and eigenfunction scale as A'"^. If a{N/2) and a{N) are 
the eigenfrequencies obtained from solutions with N/2 and 
N meshpoints, the leading-order error term therefore cancels 
in 

aj^ = ^[4a{N)-ai^N)]. (34) 

This procedure, known a s Richardson extrapolation, was 
used by Shibah ashi and Os aki (1981). It provides an esti- 
mate of the eigenfrequency that is substantially more accu- 
rate than o{N), although of course at some added computa- 
tional expense. Indeed, since the error in the representation 
|29] depends only on even powers of A^^^ the leading term of 
the error in a^i is iff{N^^). 

Even with these techniques the precision of the com- 
puted frequencies may be inadequate if the mesh used in 
stellar-evolution calculations is used also for the computa- 
tion of the oscillations. The number of meshpoints is typi- 
cally relatively modest and the distribution may not reflect 
the requirem ent to resolve properly the eigenfunct ions of 
the modes. IChristensen-Dalsgaard and Berthomieul (Il99lh 
discussed techniques to redistribute the mesh in a way that 
takes into account the asymptotic behaviour of the eigen- 
functions; a code to do so, based on four-point Lagrangian 
interpolation, is included in the ADIPLS distribution pack- 
age. On the other hand, for computing low-order modes (as 
are typically relevant for, say, 5 Scuti or j3 Cephei stars), the 
original mesh of the evolution calculation may be adequate. 

It is difficult to provide general recommendations con- 
cerning the required number of points or the need for redis- 
tribution, since this depends strongly on the types of modes 
and the properties of the stellar model. It is recommended to 

^ The variational principle is exact, formally, when the surface La- 
grangian pressure perturbation is set to zero, but not when the match to 
an isothermal atmosphere is used. 



6 



J0rgen Christensen-Dalsgaard 



carry out experiments varying the number and distribution 
of points to obtain estimates of the intrinsic precision of the 
computation (e.g.. 'Christensen-Dalsgaard an d Berthomieul 
[l991: Christensen-Dalsgaard and Mullan, 1991). In the lat- 
ter case, considering simple polytropic models, it was found 
that 4801 points yielded a relative precision substantially 
better than 10~^ for high-order p modes, when Richardson 
extrapolation was used. 

In the discussion of the frequency calculation it is im- 
portant to distinguish between precision and accuracy, the 
latter obviously referring to the extent to which the com- 
puted frequencies represent what might be considered the 
'true' frequencies of the model. In particular, the manipula- 
tions required to derive Eq.|33]and to demonstrate its varia- 
tional property depend on the equation of hydrostatic sup- 
port being satisfied. If this is not the case, as might well 
happen in an insufficiently careful stellar model calculation, 
the value determined from the variational principle may be 
quite precise, in the sense of numerically stable, but still un- 
acceptably far from the correct value. Indeed, a comparison 
between avar and CTri provides some measure of the reliabil- 
ity of the computed fr equencies (e.g. Christensen-Dalsgaard 
and Berthomieu. .l99li) . 



3 Computed quantities 



- sign 



dyi 



(35) 



Here the sum is over the zeros } in yi (excluding the cen- 
tre), and sign is the sign function, sign (z) = 1 if z > and 
sign (z) = — 1 if z < 0. For a complete model that includes 
the centre no = 1 for radial oscillations and no = for nonra- 
dial oscillations. Thus the lowest-order radial oscillation has 
order n= I. Although this is contrary to the commonly used 
convention of assigning order to the fundamental radial 
oscillation, the convention used here is in fact the more rea- 
sonable, in the sense that it ensures that n is invariant under 
a continuous variation of / from to 1 . With this definition 
n>0 for p modes, n = for f modes, and n < for g modes, 
at least in simple models. 

It has been found that this procedure has serious prob- 
lems for dipolar modes in centr ally condensed models 
(e.g., Eeei [l985: Guenth eil Il99ll: Christensen-Dalsgaard 
and Mullan, 1994). The eigenfunctions yi are shifted such 
that nodes disappear or otherwise provide spurious results 
when Eq. [35] is used to determine the mode order. A pro- 
cedure that does no t suffer from th is difficulty has recently 
been developed by iTakatal (l2006bl) : I discuss it further in 
SectionlH 

A powerful measure of the characteristics of a mode is 
provided by the normalized inertia. The code calculates this 



as 

!^^[^} + Kl+mWAr 

M[^,(/?phot)2 + Z(Z+l)^h(^phot)2] 

" 47r[yi(xphot)2+y2(xphot)V^(^ + l)] ■ ^^^^ 
(For radial modes the terms in y2 are not included.) Here 
r\ — Rx\ and = Rx^ are the distance of the innermost 
mesh point from the centre and the surface radius, respec- 
tively, and Xphot = /?phot //? = 1 is the fractional photospheric 
radius. The normalization at the photosphere is to some ex- 
tent arbitrary, of course, but reflects the fact that many radial- 
velocity observations use lines formed relatively deep in the 
atmosphere. A more common definition of the inertia is 

^mode 



E = 4kE = 



M 



(37) 



where M^ode is the so-called mode mass. 

The code has the option to output the eigenfunctions, in 
the form of {yjix")}. In addition (or instead) the displace- 
ment eigenfunctions can be output in a form indicating the 
region where the mode predominantly resides, in an energet- 
ical sense, as 



ziix) 



The programme finds th e order of the mo de a ccording to th e Z2 {x) = 
definition developed bv lScuflairel (Il974 and lOsakil Il975h . 
based on earlier work by lEckartl(ll960h . Specifically, the or- 
der is defined by = 



/ 47rr^p 
V M 

1 



1/2 



/ 4nr^p\^^^ t,r{r) 



V M 



R 



(38) 



(for radial modes only zi is found). These are defined in such 
a way that 

mz\+zl]Ax/x 

47r[yi(xphot)2+3;2(xphot)V/(/ + l)] ' 
The form provided by the z/ is also convenient, e.g., for 
computing ro tational splittings 5(U„/,„ = (!)„/„, — (0„io (e.g., 
lGoughLlT98ll) . where CO„i,„ is the frequency of a mode of ra- 
dial order n, degree I and azimuthal order m. For slow rota- 
tion the splittings are obtained from first-order perturbation 
analysis as 

5(0„i,„ = m rK„,,„{r, 0)X2(r, d)rdrdd , (40) 



characterized by kernels K„im, where in general the angular 
velocity Q depends on both r and 0. The code has built in the 
option to compute kernels for first-order rotational splitting 
in the special case where Q depends only on r. 



4 Further developments 

Several revisions of the code have been implemented in pre- 
liminary form or are under development. A substantial im- 
provement in the numerical solution of the oscillation equa- 
tions, particularly for high-order modes, is the installation of 
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a fourth-order int egratio n scheme, based on the algorithm of 
ICash and Moor3 (I1980h . This is essentially operational but 
has so far not been carefully tested. Comparisons with the re- 
sults of the variational expression and the use of Richardson 
extrapolation, of the same formal order, will be particularly 
interesting. 

As discussed by iMoya et all (l2007h the use of p' (or, as 
here, ^h) as one of the integration variables has the disad- 
vantage that the quantity A enters into the oscillation equa- 
tions. In models with a density discontinuity, such as results 
if the model has a growing convective core and diffusion is 
neglected, A has a delta-function singularity at the point of 
the discontinuity. In the ADIPLS calculations this is dealt 
with by replacing the discontinuity by a very steep and well- 
resolved slope. However, it would obviously be an advan- 
tage to avoid this problem altogether. This can be achieved 
by using instead the Lagrangian pressure perturbation 5p 
as one of the variables. Implementing this option would be 
a relatively straightforward modification to the code and is 
under consideration. 

The proper classification of dipolar modes of low or- 
der in centrally condensed models has been a long-standing 
problem in the theory of stellar pulsations, as discussed in 
Section [3] Such a scheme must provide a unique order for 
each mode, such that the order is invariant under continuous 
changes of the equilibrium model, e.g., as a result of stellar 
evolution. As a major breakthrough, Takata in a series of pa- 
pers has elucidated important properties of these modes and 
define d a new classificatio n scheme satisfying this require- 
ment dTakatal I2OO5L l2006alb '). A preliminary version of this 
scheme has been implemented and tested; however, the lat- 
est and most convenient form of the Takata classification still 
needs to be installed. 

A version of the code has been established which com- 
putes the first-order rotational splitting for a given rota- 
tion profile i2(r), in addition to setting up the correspond- 
ing kernels. This is being extended by K. Burke, Sheffield, 
to cover also sec ond-order effects of ro t ation, based on 
the formalism of iGough and ThompsonI ( Il990l) . An im- 
portant motivation for this is the integration, discussed by 
IChristensen-Dals gaardi (l2007h . of the pulsation calculation 
with the ASTEC evolution code to allow full calculation of 
oscillation frequencies for a model of specified parameters 
(mass, age, initial rotation rate, etc.) as the result of a single 
subroutine call. 
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